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1.  INTRODUCTION. 

Many  initial-boundary  value  problems  for  time-dependent  partial  differential  equa¬ 
tions  involve  fine-scale  structures  that  develop,  propagate,  decay,  and/or  disappear  as  the 
solution  evolves.  Some  examples  are  shock  waves  in  compressible  flows,  boundary  and 
shear  layers  in  viscous  flows,  and  reaction  zones  in  combustion  processes.  The  numerical 
solution  of  these  problems  is  usually  difficult  because  the  nature,  location,  and  duration  of 
the  structures  are  often  not  known  in  advance.  Thus,  conventional  numerical  approaches 
that  calculate  solutions  on  a  prescribed  (typically  uniform)  mesh  often  fail  to  adequately 
resolve  the  fine-scale  phenomena,  have  excessive  computational  costs,  or  produce  incorrect 
results.  Adaptive  procedures  that  evolve  with  the  solution  offer  a  robust,  reliable,  and 
efficient  alternative.  Such  techniques  have  been  the  subject  of  a  great  deal  of  recent  atten¬ 
tion  (cf.  Babuska  et  al.  [7,  9])*  and  are  generally  capable  of  introducing  finer  meshes  in 
regions  where  greater  resolution  is  needed  [1,  2,  3,  6,  8,  10,  15,  16],  moving  meshes  in 
order  to  follow  isolated  dynamic  phenomena  [1,  2,  5,  21,  23,  24,  25,  30],  or  changing  the 
order  of  methods  in  specific  regions  of  the  problem  domain  [18,  22].  The  utility  of  such 
adaptive  techniques  is  greatly  enhanced  when  they  are  capable  of  providing  an  estimate  of 
the  accuracy  of  the  computed  solution.  Local  error  estimates  are  often  used  as  refinement 
indicators  and  to  produce  solutions  that  satisfy  either  local  or  global  accuracy 
specifications  [1,  2,  3,  6,  8,  10,  15,  16].  Successful  error  estimates  have  been  obtained 


•  References  are  listed  at  the  end  of  this  report. 


using  h-refinement  [6,  15,  16],  where  the  difference  between  solutions  on  different  meshes 


is  used  to  estimate  the  error,  and  p-iefinement  [1,  2,  3,  8,  16,  22]  where  the  differences 
between  methods  of  different  orders  are  used  to  estimate  the  error. 


We  discuss  an  adaptive  procedure  that  combines  mesh  movement  and  local 
refinement  for  m-dimensional  vector  systems  of  partial  differential  equations  having  the 
form 


u,  +  f(a:,y,r,u,u^,Uj,)  =  +  [D^(:c,y,r,u)Uj,]y, 

for  f  >  0,  ix,y)  s  Q, 


(la) 


with  initial  conditions 


u(x ,y  ,0)  =  ,y ),  for  (;c  ,y )  6  Q  vj  8Q,  (lb) 

and  appropriate  well-posed  boundary  conditions  on  the  boundary  dCl  of  a  rectangular 
region  Q. 


We  suppose  that  a  numerical  method  is  available  for  calculating  approximate  solu¬ 
tions  and  error  indicators  of  Eq.  (1)  at  each  node  of  a  moving  mesh  of  quadrilateral  cells. 
Any  appropriate  numerical  method  is  applicable  and  the  error  indicator  can  either  be  an 
estimate  of  the  local  discretization  error  or  another  function  (e.g.,  an  estimate  of  the  solu¬ 
tion  gradient  or  curvature)  that  is  large  where  additional  resolution  is  needed  and  small 
where  less  resolution  is  desired.  Our  adaptive  algorithm  consists  of  three  main  parts:  (a) 
movement  of  a  coarse  base  mesh,  (b)  local  refinement  of  the  base  mesh  in  regions  where 
resolution  is  inadequate,  and  (c)  creation  and  regeneration  of  the  base  mesh  when  it 
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becomes  overly  distorted.  Our  experience  (cf.  Section  III)  and  that  of  others  [25]  indicates 
that  mesh  motion  can  substantially  reduce  errors  for  a  very  modest  computational  cost. 
Mesh  motion  alone,  however,  cannot  produce  a  solution  that  will  satisfy  a  prescribed  error 
tolerance  in  all  situations.  For  this  reason,  we  have  combined  mesh  motion  with  local 
mesh  refinement,  and  recursively  solve  local  problems  in  regions  where  error  tolerances 
are  not  satisfied.  The  local  solution  scheme  successively  reduces  the  domain  size  and, 
thus,  further  reduces  the  cost  of  the  computation.  Some  problems,  e.g.,  those  with  severe 
material  deformations,  can  result  in  tangling  and  distortion  of  the  moving  base  mesh. 
Therefore,  we  have  created  a  procedure  that  automatically  generates  a  new  base  mesh 
whenever  the  old  one  is  unsuitable. 

The  adaptive  procedures  described  in  this  report  combine  our  earlier  work  on  mesh 
moving  techniques  [5]  and  local  refinement  procedures  [6].  The  inclusion  of  a  static  mesh 
regeneration  scheme  adds  greater  reliability  and  efficiency  to  these  methods.  The  three 
components  of  our  adaptive  algorithm  are  described  in  Section  11;  however,  frequent  refer¬ 
ences  are  made  to  our  previous  investigations  [5,  6].  A  computer  code  based  on  the  adap¬ 
tive  algorithm  of  Section  II  has  been  combined  with  a  MacCormack  finite  difference 
scheme  and  an  error  indicator  based  on  Richardson  extrapolation.  It  has  been  used  to 
solve  a  sequence  of  hyperbolic  problems  (i.e.,  problems  having  the  form  of  Eq.  (1)  with 
D'  =  :=  0)  and  our  findings  on  three  examples,  where  we  have  attempted  to  appraise 

the  relative  costs  and  benefits  of  the  mesh  moving  and  local  refinement  portions  of  our 


adaptive  algorithm,  are  reported  in  Section  III.  We  have  also  compared  solutions  obtained 
by  adaptive  techniques  to  those  obtained  using  stationary  uniform  meshes.  In  all  three 
examples,  solutions  obtained  by  adaptive  techniques  cost  less  than  solutions  obtained  on 
stationary  uniform  meshes  having  approximately  the  same  accuracy.  The  mesh  moving 
technique  added  approximately  ten  percent  to  the  computational  time  of  the  adaptive  algo¬ 
rithm  and  greatly  improved  the  results.  Most  of  the  computational  time  was  devoted  to 
calculating  the  solution  and  error  indicators,  and  not  to  the  overhead  induced  by  the 
refinement  procedure.  Although  we  are  greatly  encouraged  by  our  results,  our  adaptive 
procedures  are  far  from  complete.  Some  possible  improvements  and  future  considerations 
are  discussed  in  Section  FV. 

II.  ALGORITHM  DESCRIPTION. 

A  top-level  description  of  our  adaptive  procedure  is  presented  in  Figure  1  in  a 
pseudo-PASCAL  language.  This  procedure  is  called  adaptive_PDE_solver  and  it 
integrates  a  system  of  partial  differential  equations  fi’om  time  tinit  to  tfinal  and  attempts 
to  keep  the  local  error  indicators  below  a  tolerance  of  tol.  The  base  level  time  step  Ar  is 
initially  specified,  but  may  be  changed,  as  needed,  during  the  integration. 

The  rectangular  domain  Q  is  initially  discretized  into  a  coarse  moving  spatial  grid  of 
M  xN  quadrilateral  cells.  An  initial  base  mesh  is  generated  from  this  mesh  by  increasing 
the  values  of  M  and  N ,  as  necessary,  and  moving  the  mesh  so  that  it  is  concentrated  in 


procedure  a(iaptive_PDE_soIver(ri/ur,  At,  tfinal,  tol:  real;  M,N:  integer); 


begin 

Generate  an  initial  base  mesh; 
t  :=  tinif, 

while  t  <  tfinal  do 
begin 

Move  the  base  mesh  for  the  time  step  t  to  r  +  Ar; 
local_refine(0,  t.  At,  fol); 
t  :=  t  +  At; 

Select  an  approprir  .e  At; 

if  base  mesh  is  too  distoned  then  regenerate  a  base  mesh 
end 

end  {  adaptive_PDE_solver  ); 


Figure  1.  Pseudo-PASCAL  description  of  an  adaptive  procedure  to  solve  the 
partial  differential  system  in  Eq.  (1)  from  t  =  tinit  to  tfinal  to  within  a  toler¬ 
ance  of  tol . 


regions  where  error  indicators  are  large  (cf.  Section  11.3).  The  base  mesh  is  moved  for 
each  base  time  step  Ar  (cf.  Amey  and  Flaherty  [5]  and  Section  11. 1)  and  the  partial 
differential  system  in  Eq.  (1)  is  solved  on  this  mesh  for  a  base  time  step.  This  is  followed 
by  a  recursive  local  mesh  refinement  in  regions  where  error  indicators  are  larger  than  tol . 
The  local  refinement  procedure  localjrefine  was  described  in  Amey  and  Flaherty  [6]  and 
its  major  features  are  summarized  in  Section  II.2.  The  integration  for  each  base-mesh 
time  step  is  concluded  by  the  selection  of  a  new  value  of  At  for  the  subsequent  time  step 


and  the  generation  of  new  base  mesh  (cf.  Section  n.3),  if  necessary. 

The  mesh  moving,  local  refinement,  and  mesh  regeneration  algorithms  are  uncoupled 
from  each  other  as  well  as  from  the  procedures  used  to  solve  the  partial  differential  system 
and  calculate  local  error  indicators.  This  reduces  computational  costs  and  provides  a  great 
deal  of  flexibility.  Thus,  individual  modules  can  easily  be  replaced,  omitted,  or  combined 
with  other  software. 

ILl.  Mesh  Moving  Algorithm.  Mesh  moving  strategies  should  produce  a  smooth  mesh 
where  the  sizes  of  neighboring  computational  cells  vary  slowly  and  cell  angles  differ  only 
by  modest  amounts  from  right  angles.  It  is,  of  course,  essential  for  the  nodes  of  the  mesh 
to  remain  within  Q  and  for  cells  not  to  overlap.  Meshes  that  violate  these  conditions  can 
produce  large  discretization  errors  that  overwhelm  the  positive  effects  of  mesh  moving. 
Our  mesh  moving  procedure  is  based  on  an  intuitive  approach  rather  than  more  analytic 
error  equidistribution  (cf.,  e.g.,  Coyle  et  al.  [19]  or  Dwyer  [23])  and  variational  approaches 
(cf.  ^rackbill  and  Saltzman  [17]).  The  essential  idea  is  to  move  the  mesh  so  that  it  may 
follow  isolated  nonuniformities,  such  as  wave  fronts,  shock  layers,  and  reaction  zones. 
This  generally  reduces  dispersive  errors  and  allows  the  use  of  larger  time  steps  while 
maintaining  accuracy  and  stability. 

At  each  base  time,  we  scan  the  M  xA^  base  mesh  of  quadrilateral  cells  and  locate 
"significant  error  nodes"  as  those  having  error  indicators  greater  than  twice  the  mean  nodal 
error  indicator  and  also  greater  than  ten  percent  of  tol.  This  empirical  strategy  avoids 
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having  the  mesh  respond  to  fluctuations  when  error  indicators  are  too  small,  but  is  sensi¬ 
tive  enough  to  avoid  missing  dynamic  phenomena  associated  with  large  error  indicators. 
If  there  are  no  significant  error  nodes,  computation  proceeds  on  a  stationary  mesh.  The 
nearest  neighbor  clustering  algorithm  of  Berger  and  Oliger  [15]  is  then  used  to  gather  the 
significant  error  nodes  into  clusters.  In  this  iterative  algorithm,  a  cluster  is  first  defined  to 
consist  of  one  arbitrary  significant  error  node.  Other  significant  error  nodes  are  added  to 
the  cluster  if  they  are  within  a  specified  minimum  intercluster  distance  from  the  nearest 
node  in  the  cluster.  We  take  the  minimum  intercluster  distance  to  be  the  length  of  a  cell 
diagonal.  New  clusters  are  established  for  nodes  that  do  not  belong  to  any  existing  clus¬ 
ter.  Clusters  are  united  when  a  node  is  determined  to  belong  to  more  than  one  of  them. 
Upon  completion  of  the  algorithm,  (a)  nodes  in  different  clusters  will  be  separated  by  at 
least  the  minimum  intercluster  distance,  and  (b)  no  node  in  a  cluster  with  more  than  one 
node  will  be  further  than  the  minimum  intercluster  distance  from  its  nearest  neighbor  in 
the  cluster. 


Following  Berger  and  Oliger  [15],  we  generate  near  minimum  area  rectangles  that 
contain  each  cluster.  The  principal  axes  of  each  rectangle  are  the  major  and  minor  axes  of 
an  enclosed  ellipse  having  the  same  first  and  second  moments  as  the  nodes  in  the  cluster. 
Thus,  if  U,  )  are  the  coordinates  of  a  node  and  ,y^ )  are  the  mean  coordinates  of  all 
nodes  in  the  cluster,  then  the  axes  of  the  rectangle  are  in  the  directions  of  the  eigenvectors 
of  the  symmetric  (2x2)  matrix 


E  -  x^)  £  (Xiyi  -  x„y„ ) 

E  (Xiyi  -  x^y^)  ^  (y,-2  -  y^) 


(2) 


L  J 

where  the  summations  range  over  all  nodes  in  the  cluster. 

For  many  problems,  there  may  be  too  small  a  percentage  of  significant  error  nodes 
within  a  cluster.  In  order  to  reduce  this  inefficiency  and  provide  some  alignment  with, 
e.g.,  curved  wave  fronts,  the  rectangles  are  checked  for  efficiency  by  determining  the  per¬ 
centage  of  significant  error  nodes  in  each  rectangle.  If  a  fifty-percent  efficiency  is  not 
achieved,  the  rectangle  is  iteratively  bisected  in  the  direction  of  its  major  axis  until  all 
clusters  have  at  least  a  fifty-percent  efficiency. 


We  determine  node  movement  from  the  velocity  of  propagation,  the  orientation,  and 
the  size  of  error  clusters.  We  assume  that  nodes  in  the  same  cluster  have  related  solution 
characteristics,  so  that  we  can  determine  individual  node  movement  from  the  propagation 
of  the  center  of  the  error  cluster.  Each  cluster  moves  according  to  the  differential  equa¬ 
tion 

Tm  +  =  0.  (3) 

where  r„(t)  -  [x„({),y„(t)]^  is  the  position  of  the  center  of  an  error  cluster  ^d 
(')  ;=  d(  )/dt.  The  choice  of  the  parameter  X  can  be  critical  in  certain  situations.  If  k  is 
selected  too  large,  the  system  in  Eq.  (3)  will  be  stiff  and  computationally  expensive.  On 
the  other  hand,  if  X  is  too  small,  the  mesh  can  oscillate  from  time  step-to-time  step. 
Coyle  et  al.  [19]  and  Adjerid  and  Flaherty  [2]  suggested  some  adaptive  procedures  for 


choosing  X;  however,  we  found  no  appreciable  differences  in  results  or  computation  times 
when  X  varied  significantly.  The  examples  of  Section  HI  were  calculated  with  X  =  1. 

We  solve  Eq.  (3)  for  each  base  time  step  and  each  cluster  using  an  explicit  numerical 
method.  The  center  of  an  error  cluster  is  moved  a  distance  Ar^  =  r^_(r+Ar)  -  T„it)  at 
the  base  time  f.  Let  and  denote  the  projections  of  Ar^  onto  the  major  and 
minor  axes  of  the  rectangular  cluster.  We  use  the  one-dimensional  piecewise  linear  func¬ 
tion 


Ar^(3/2+x,7w,  ),  if -3w,  /2  ^  a:,-  ^ -w,-/2 
Ar^. ,  if  -w;  /2  <  Xi  <  wi  12 

Ar^i3l2-XilWi),  if  w,/2  ^  Xi  ^  3w',-/2 


(4) 


[  0,  otherwise 

to  move  the  nodes  of  the  mesh  along  the  two  principal  axial  directions  of  the  error  clus¬ 
ters.  The  cluster  referred  to  in  Eq.  (4)  has  dimensions  >viXW2,  and  (xi,X2)  are  local 
Cartesian  coordinates  of  a  node  in  the  principal  directions  of  the  cluster  relative  to  its 
center.  For  i  =  1,  nodes  in  the  range  of  the  cluster  (-3w  i/2  $  x  i  ^  3w  i/2, 
-W2/2  <  X2  ^  >^2/2)  are  moved  a  distance  dj  This  situation  is  shown  in  Figure  2. 


In  order  to  maintain  smooth  mesh  motion  throughout  the  domain,  nodes  outside  the 
range  of  a  cluster  move  in  a  distance 


^i,oiusuU  ~  inside  ~  (2r/Z))],  i  =  1,  2, 


(5) 


s 


domain  Q 


Figure  2.  A  rectangular  wiXW2  error  cluster.  Nodes  within  the  range  of  the 
cluster,  Swi  X  W2,  are  moved  a  distance  di  i,^  in  the  jcj  principal  direction  ac¬ 
cording  to  Eq.  (4).  Nodes  outside  the  range  of  the  cluster  are  moved  a  distance 
d  1  outside  in  the  X I  direction  according  to  Eq.  (5).  The  distance  z  is  the  shortest 
distance  to  the  range  of  the  cluster. 


where  z  is  the  shonest  distance  to  the  range  of  the  cluster  (cf.  Figure  2)  and  D  is  the 
diagonal  of  Q.  For  each  cluster,  the  mesh  is  moved  in  the  direction  of  the  major  axis 
(i  =  1)  using  Eqs.  (4)  and  (5).  This  is  followed  by  a  similar  procedure  in  the  direction  of 
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the  minor  axis  (i  =  2).  The  distances  and  di  outside  ^  reduced  near  3Q  in  order 

to  prevent  nodes  from  leaving  fl.  In  particular,  we  recalculate  di^  as  dij[min(\,  bic)}, 
i  =  \,  2,  j  =  inside,  outside,  where  b  is  the  distance  of  the  node  to  the  boundary  and  c  is 
twice  the  length  of  a  cell  diagonal  on  a  uniform  mesh  having  the  same  number  of  cells  as 
the  moving  mesh.  Nodes  on  domain  boundaries,  except  comer  nodes,  which  are  not  ^ 
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moved,  are  restrained  to  move  along  the  boundary.  Finally,  the  mesh  moving  algorithm  is 


not  restricted  to  the  functions  given  by  Eqs.  (4)  and  (5),  and  several  other  choices  are  pos¬ 


sible. 


II.2.  Local  Refinement  Algorithm,  As  shown  in  Figure  1,  the  local  refinement  procedure 


is  invoked  after  the  base  mesh  has  been  moved  for  a  base  time  step.  Our  refinement  stra¬ 


tegy  consists  of  first  calculating  a  preliminary  solution  on  the  base  mesh  for  a  base  time 


step.  An  error  indicator  is  used  to  locate  regions  where  greater  resolution  is  needed. 


Finer  grids  are  adaptively  created  in  these  high-error  regions  by  locally  bisecting  the  time 


step  and  the  sides  of  the  quadrilateral  cells  of  the  base  grid  and  the  solution  and  error  indi¬ 


cators  are  computed  on  the  finer  grids.  The  refinement  scheme  is  recursive;  thus,  fine 


subgrids  may  be  refined  by  adaptively  creating  even  finer  subgrids.  This  relationship  leads 


naturally  to  a  tree-data  structure.  Information  regarding  the  geometry,  solution,  and  error 


indicators  of  the  base  grid  is  stored  as  the  root  node  or  level  0  of  the  tree.  Subgrids  of  the 


base  grid  are  offsprings  of  the  root  node  and  arc  stored  as  level  1  of  the  tree.  The  struc¬ 


ture  continues,  with  a  grid  at  level  /  having  a  parent  coarser  grid  at  level  /  -  1  and  any 


EV-V 


finer  offspring  grids  at  level  /  +  1.  Grids  at  level  /  of  the  tree  are  given  an  arbitrary  ord¬ 
ering  and  we  denote  them  as  Gij,  j  =  1,  2, ....  N/,  where  A//  is  the  number  of  grids  at 
level  /.  Our  refinement  procedures  permit  grids  at  the  same  level  of  a  two-dimensional 
problem  to  intersect  and  overlap;  however  offspring  grids  must  be  properly  nested  within 
the  boundaries  of  their  parent  grid.  A  one-dimensional  grid  with  its  appropriate  4ree  struc¬ 
ture  for  a  base  rime  step  is  shown  in  Figure  3. 

A  top-level  pseudo- PAS  CAL  description  of  a  recursive  local  refinement  algorithm 
that  solves  systems  of  the  form  in  Eq.  (1)  on  the  tree  of  grids  described  above  is  presented 
in  Figure  4.  The  procedure  local_refine  integrates  partial  differential  equations  on  the 
grids  Gi^j,  j  =  1,  2,  Ni,  at  level  /  of  the  tree  from  time  tinit  to  timt  +  A/  and  attempts 
to  satisfy  a  prescribed  local  error  tolerance  tol.  For  each  grid  at  level  /,  a  solution  and 
error  indicators  are  calculated  at  rime  tinit  +  Ar.  Additional  finer  grids  are  introduced  in 
regions  where  the  error  indicators  exceed  the  prescribed  tolerance  tol  and  the  differential 
system  is  solved  again  on  the  finer  grids  using  two  rime  steps  of  duration  Ar/2  and  a  toler¬ 
ance  of  tolH.  Observe  that  the  solution,  error  indicators,  and  refined  subgrids  are  calcu¬ 
lated  for  all  grids  at  level  /  before  calculating  any  solutions  at  level  /  +  1.  Implicit  in 
local_refine  are  the  assumptions  that  a  solution  can  be  computed  on  any  grid  and  that 
refinement  terminates.  If  either  of  these  assumptions  are  violated,  the  procedure  ter¬ 
minates  in  failure. 

Our  technique  for  introducing  finer  subgrids  consists  of  four  steps:  (a)  an  initial  scan 


Figure  3.  Coarse  and  refined  grids  (top)  and  their  tree  representation  (bottom) 
for  a  one-dimensional  example. 


of  each  level  /  grid  to  locate  "untolerable-error"  nodes  as  those  where  the  error  indicator 
exceeds  the  prescribed  tolerance  tol,  (b)  clustering  any  untolerable  nodes  into  rectangular 
regions,  (c)  buffering  the  clustered  regions  in  order  to  reduce  problems  associated  with 
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procedure  local_refme(  /:  integer;  tinit,  tol:  real ); 
begin 

for  j  :=  1  to  N[l]  do 
begin 

Integrate  the  partial  differential  system  from  tinit  to  Unit  +  A/ 
on  grid  G[IJ]; 

Calculate  error  indicators  at  tinit  +  A/  at  all  nodes  of  grid 

GUJV, 

if  any  error  indicators  >  tol  then  introduce  level  /  -i-  1  subgrids 
ofG[/j] 
end  (  for  }; 

if  any  error  indicators  >  tol  then 
begin 

local_refine(/  +  1,  tinit,  At/2,  tol/2); 

local_refine(/  +  1,  tinit  +  Ar/2,  At/2,  tol/2) 
end 

end  {  local_refine  }; 


Figure  4.  Pseudo-PASCAL  description  of  a  recursive  local  refinement  pro¬ 
cedure  to  find  a  solution  of  the  partial  differential  system  in  Eq.  (1)  on  all  grids 
at  level  /  of  the  tree. 


prescribing  initial  and  boundary  conditions  at  coarse/fine  grid  interfaces,  and  (d)  cellularly 
refining  the  level  /  meshes  and  time  step  within  the  buffered  clusters.  Of  course,  if  there 
are  no  untolerable-error  nodes,  the  solution  is  acceptable  and  further  refinement  is 


unnecessary. 


The  same  clustering  algorithm  of  Berger  and  Oliger  [15]  that  was  used  to  move  the 


base  mesh  is  also  used  to  group  untolerable-error  nodes  for  refinement.  Each  rectangular 
error  cluster  is  enlarged  by  increasing  its  major  and  minor  axes  by  twice  the  size  of  the 
average  ceil  edge  within  the  cluster.  The  region  between  the  enlarged  and  original  error 
clusters  provides  a  buffer  so  that  artificial  internal  boundary  conditions  (that  are  discussed 
below)  will  be  prescribed  at  low-error  nodes  as  far  as  possible  and  fine-grid  errors  will  not 
propagate  through  the  buffer  in  a  time  step. 

Refined  subgrids  are  created  by  bisecting  the  time  step  and  edges  of  each  cell  of  the 
parent  mesh  that  intersects  the  buffered  rectangular  error  clusters.  Coarse  mesh  motion  is 
maintained  on  the  refined  grids  so  that  after  two  time  steps  of  size  Ar/2,  cells  of  the 
refined  grids  will  be  properly  nested  within  those  of  their  parent  grid.  Additional  details 
of  the  refinement  algorithm  and  data  structures  are  presented  in  Amey  and  Flaherty  [6]. 

Artificial  initial  and  boundary  data  must  be  determined  from  solutions  on  other  grids 
in  order  to  calculate  the  solution  and  error  indicators  on  refined  subgrids.  Furthermore, 
solutions  on  finer  grids  are  used  to  replace  those  on  coarser  grids  at  common  nodal  loca¬ 
tions. 

Initial  data  for  a  subgrid  are  calculated  directly  from  the  initial  function  uVj)  at 
r  =  0.  For  t  >  0,  initial  data  are  obtained  by  interpolation  using  the  solution  at  the  same 
time  on  the  finest  available  mesh.  In  order  to  provide  data  for  this  interpolation,  we  save 
all  solution  values  on  previous  subgrids  until  they  are  no  longer  needed  due  to 
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advancement  in  time  of  an  acceptable  solution.  Bilinear  functions  using  the  solution 
values  at  the  four  vertices  of  the  finest  existing  cell  are  used  to  obtain  the  solution  at  the 
nodes  of  cells  of  the  refined  mesh.  Further  analysis  is  needed  regarding  the  effects  on 
accuracy  and  stability  and  the  proper  order  of  this  interpolation.  Bieterman,  Raherty,  and 
Moore  [16]  give  an  example  where  the  fine-scale  structure  of  a  solution  was  lost  by  inter¬ 
polation  from  too  coarse  a  mesh. 

In  a  similar  manner,  boundary  data  for  refined  meshes  are  calculated  directly  from 
the  prescribed  boundary  conditions  on  portions  of  subgrids  that  intersect  dCl.  Dirichlet 
boundary  data  are  prescribed  on  the  edges  of  subgrids  that  are  in  the  interior  of  by 
interpolating  the  solution  from  coarser  meshes.  Bilinear  functions  using  the  solution 
values  at  the  four  vertices  of  the  adjacent  face  of  the  finest  existing  space-time  cell  are 
used  to  obtain  solution  values  for  the  nodes  of  refined  cells. 

Acceptable  fine-mesh  solutions  are  used  to  replace  solutions  at  the  nodes  of  coarser 
grids  that  lie  within  the  untolerable-error  portions  of  clusters.  Solutions  at  low-error  nodes 
in  the  buffer  zones  of  clusters  are  not  replaced  in  order  to  avoid  possible  contamination  of 
accurate  solutions.  When  fine  grids  overlap  each  other  in  an  untolerable-error  region,  the 
average  value  of  the  solutions  at  common  fine-grid  nodes  is  used  to  replace  the  appropri¬ 
ate  coarse  grid  solution.  Boundary  effects  do  not  propagate  through  a  sufficiently  large 
buffer  and,  thus,  have  no  effect  on  the  solution  within  the  untolerable-error  region  of  a 
cluster  when  an  explicit  numerical  scheme  is  used  for  the  integration.  Greater  care  is 


needed  when  implicit  integration  methods  are  used,  since  artificial  boundary  conditions 
can  affect  the  accuracy,  convergence,  and  stability  of  the  solution  at  all  nodes  in  the  clus¬ 
ter  regardless  of  the  size  of  the  buffer. 

Stability  and  conservation  of,  e.g.,  fluxes  at  interfaces  between  coarse  and  fine 
meshes  must  be  investigated  further,  particularly  in  two  dimensions.  For  one-dimensional 
problems,  Berger  and  Oliger  [15]  showed  that  linear  interpolation  of  solutions  from  a 
coarse  to  a  fine  mesh  produced  no  instabilities  in  the  the  Lax-Wendroff  scheme.  Berger 
[14]  also  discussed  conservation  at  mesh  interfaces  and  proposed  explicit  enforcement  of 
conserved  quantities  at  coarse/fine  mesh  boundaries.  Rai  [29]  presented  some  finite 
difference  schemes  that  maintained  conservation  at  grid  interfaces  for  two-dimensional 
compressible  flow  problems. 

IL3.  Initial  Mesh  Construction  and  Regeneration.  The  efficiency  of  our  adaptive  mesh 
moving  and  refinement  strategies  is  dependent  on  our  ability  to  generate  a  suitable  initial 
mesh  and  to  regenerate  a  new  base  mesh  should  it  become  severely  distorted  at  later 
times.  The  proper  base  mesh  can  reduce  the  need  for  refinement  and,  thus,  increase 
efficiency. 

The  two  essential  elements  of  a  mesh  generation  or  regeneration  procedure  are  the 
determination  of  the  number  of  nodes  and  their  optimal  location.  A  base  mesh  having  too 
few  nodes  will  result  in  excessive  refinement  while  one  having  too  many  nodes  will 
reduce  efficiency.  Many  mesh  generation  procedures  have  been  developed  (cf.,  e.g., 


Thompson  [31]  or  Brackbill  and  Saltzman  [17]);  however,  the  best  one  to  use  in  conjunc¬ 
tion  with  an  adaptive  procedure  is  still  far  from  being  established.  Our  current  approach 
to  mesh  generation  is  to  use  the  error  indicators  computed  by  a  trial  solution  to  determine 
an  initial  mesh  that  approximately  equidistributes  the  error  indicators. 

To  begin,  we  create  a  uniform  M  xN  rectangular  mesh  using  prescribed  values  of  M 
and  N  that  reflect  the  coarsest  mesh  that  should  be  used  to  calculate  a  solution.  We  solve 
the  system  in  Eq.  (1)  for  a  base  time  step  Ar  on  the  uniform  stationary  base  mesh  and 
compute  the  solution  and  error  indicators.  Local  mesh  refinement  is  performed  as 
described  in  Section  n.2  until  the  prescribed  tolerance  is  attained.  We  use  this  solution  to 
determine  the  number  of  nodes  AT  in  a  new  base  mesh  as 

a:  =M1V  +  £(3/4)'a:,,  (6a) 

1=1 

where  ATj  is  the  number  of  nodes  introduced  at  level  /  and  n  is  the  total  number  of  levels 
in  the  tree.  Having  computed  K ,  we  calculate  the  dimensions  of  a  new  M  xN  mesh  as 

M  =  <kWN,  N  =  WTJm.  (6b) 

The  bars  have  been  omitted  on  M  and  N  in  the  algorithms  displayed  in  Figures  1  and  4 
and  in  all  further  discussions. 

Node  placement  for  the  new  base  mesh  is  accomplished  by  locating  all  nodes  of  the 


original  base  mesh  having  error  indicators  that  are  greater  than  twice  the  mean  error  indi- 


cator.  These  ncxles  are  then  grouped  into  rectangular  clusters  using  the  clustering  algo¬ 
rithm  of  Section  11,1.  A  uniform  base  mesh  is  generated  when  there  are  no  nodes  having 
error  indicators  that  are  greater  than  twice  the  mean  error  indicator. 

Nodes  are  moved  towards  the  center  of  the  nearest  error  cluster  unless  they  are 
within  a  two-cell  diagonal  range  of  two  or  more  error  clusters.  In  the  former  case,  a  node 
is  moved  four-tenths  of  its  distance  to  the  center  of  the  nearest  cluster  unless  this  distance 
is  greater  than  12.5  times  the  average  cell  diagonal,  in  which  case  it  is  moved  five  times 
the  a'/erage  cell  diagonal.  Nodes  that  are  within  a  two-cell  diagonal  range  of  two  or  more 
clusters  are  moved  by  four- tenths  of  a  weighted  average  of  the  distances  to  centers  of  the 
involved  clusters.  Nodes  on  9Q  remain  on  9Q.  Nodes  near  the  boundary  move  a  reduced 
distance  in  order  to  prevent  the  formation  of  large  elements.  When  an  error  cluster  inter¬ 
sects  opposite  boundaries  of  nodes  are  not  moved  in  the  direction  of  the  major  axis  of 
the  cluster.  This  construction  generates  a  base  mesh  that  depends  on  the  solution  of  the 
panial  differential  system  as  well  as  its  initial  condition. 

The  base  mesh  can  become  severely  distorted  for  some  problems  (cf.  Amey  and 
Flaherty  [5])  and  we  would  like  a  capability  for  generating  a  new  base  mesh  whenever 
this  happens.  Since  the  new  mesh  is  created  at  a  specific  time,  rather  than  by  mesh 
motion,  we  refer  to  this  process  as  static  mesh  regeneration.  Our  static  mesh  regeneration 
procedure  consists  of  three  steps:  (a)  determining  that  there  is  a  need  for  a  new  base  mesh, 
(b)  creating  the  new  base  mesh,  and  (c)  interpolating  the  solution  from  the  old  to  the  new 


base  mesn. 


A  mesh  is  regenerated  when  any  interior  angle  of  a  cell  is  less  than  50  or  greater 
than  130  degrees,  the  aspect  ratio  of  any  ceil  is  greater  than  15,  or  the  mesh  ratio  of  adja¬ 
cent  cells  exceeds  5  or  is  less  than  1/5.  In  the  present  context,  the  aspect  ratio  is  defined 
as  the  average  length  divided  by  the  average  width  of  a  cell  and  the  mesh  ratios  are 
defined  as  the  ratio  of  the  lengths  and  widths  of  adjacent  cell  sides. 


A  new  base  mesh,  having  the  same  number  of  nodes  as  the  old  one,  is  generated 
using  the  procedure  described  above  for  creating  an  initial  base  mesh.  The  error  clusters 
for  the  existing  mesh  are  used  to  generate  the  new  base  mesh,  so  that  new  clusters  do  not 
have  to  be  computed.  This  process  appears  to  reduce  angle  deviations  from  ninety 
degrees,  control  aspect  ratios,  and  mollify  adjacent  mesh  ratios. 


Once  a  new  base  mesh  has  been  constructed,  the  solution  on  the  old  one  is  interpo¬ 
lated  to  the  new  one  by  using  bilinear  interpolation  with  respect  to  the  cells  of  the  old 
base  mesh.  The  order  and  nature  of  the  interpolation  needs  further  investigation  and  we 
are  studying  methods  that  conserve,  e.g.,  fluxes  (cf.  Berger  [14]  or  Rai  [29]). 


rn.  COMPUTATIONAL  EXAMPLES. 


In  order  to  demonstrate  the  capabilities  of  the  adaptive  procedure  described  in  Sec¬ 
tion  II,  we  applied  it  to  three  hyperbolic  systems.  We  used  a  two-step  MacCormack  finite 
difference  method  (cf.  Amey  and  Flaherty  [5],  Hindman  [26],  or  MacCormack  [27])  to 


integrate  the  partial  differential  equations  and  Richardson’s  extrapolation  (cf.  Amey  [4]  or 
Berger  and  Oliger  [15])  to  indicate  local  errors.  Base  mesh  geometry  was  prescribed  as 
indicated  in  each  example.  If  the  base  mesh  time  step  failed  to  satisfy  the  Courant, 
Friedrichs,  Lewy  theorem,  it  was  automatically  reduced  to  the  maximum  allowed  by  the 
Courant  condition  (cf.  Amey  [4]  and  Amey  and  Flaherty  [6]).  This  procedure  should  also 
satisfy  the  Comant  condition  on  all  subgrids  when  the  characteristic  speeds  vary  slowly. 

Numerical  results  obtained  on  uniform  stationary  grids  are  compared  with  those 
obtained  by  adaptive  strategies  that  use  (a)  mesh  moving  only,  (b)  local  refinement  only, 
and  (c)  the  combination  of  mesh  moving  and  refinement  discussed  in  Section  H.  The 
examples  are  designed  to  determine  the  relative  cost,  accuracy,  and  efficiency  of  our  adap¬ 
tive  algorithm  and  each  of  its  components.  Accuracy  is  appraised  by  computing  the 
difference  e  between  the  exact  and  numerical  solutions  of  a  problem  in  either  the  max¬ 
imum  or  L 1  norms,  i.e.,  by  computing  either 


l|e(  /,t)IL  :=  max  max  I  e .  (x,- ,y,- ,r )  I . 

l<i</:  \<jim  * 


(7a) 


or 


l|e(-.%f)|li=JJP  I:  \ej\dxdy,  (7b) 

£1  >=i 

respectively.  Here,  K  is  the  number  of  nodes  in  the  mesh  at  time  t  and  P  is  a  piecewise 
constant  interpolation  operator  with  respect  to  the  cells  of  the  base  mesh  that,  on  each  cell. 


has  the  average  value  of  the  errors  at  the  vertices  of  the  ceU.  We  use  either  the  total  CPU 


time  or  the  maximum  number  of  nodes  used  in  a  base  time  step  as  measures  of  the  com¬ 


putational  complexity  of  a  procedure.  All  calculations  were  performed  in  double  precision 
arithmetic  on  an  IBM  3081/D  computer  at  the  Rensselaer  Polytechnic  Institute. 

Solutions  are  displayed  by  drawing  either  level  lines  or  wire-frame  perspective  rendi¬ 
tions.  Meshes  are  displayed  by  showing  the  complete  two-dimensional  spatial  discretiza¬ 
tion  at  specified  times  with  finer  subgrids  overlaying  coarser  ones.  This  portrayal  does  not 
show  the  reduced  time  steps  that  are  used  for  the  subgrid  calculations.  The  broken-line 
rectangles  in  the  figures  indicate  the  error  cluster(s)  that  are  used  to  move  the  base  mesh. 

Example  1.  Consider  the  linear  initial- boundary  value  problem  that  was  proposed  as 
a  test  problem  by  McRae  et  aJ.  [28]: 


M,  -  yUjf  +  xUy  =0,  t  >  0,  (x  ,y )  6  Q,  (8a) 

Jo,  if  (x-1/2)2+ 1.5y2  >  1/16 

w(x,y,0)  1  -  16((x-1/2)^  +  1.5y^),  otherwise, 

{x,y)en[^dCl,  (8b) 

and 

u{x,y,t)  =  0,  t  >0,  {x,y)  e  dQ,  (8c) 

where  Q  :=  { (x,y)  I  -1.2  <x,y  <  12]. 

The  exact  solution  of  Eq.  (8)  is  an  elliptical  cone  that  rotates  about  the  origin  in  the 


counterclockwise  direction  witli  period  2k.  It  can  be  written  in  the  form 


where 


0,  if  C  <  0 

ifcso, 

» 


(9a) 


C  =  1  -  16[(j:cosr  +ysinr  -  1/2)^+  1.5(ycosr  -xsin/)^].  (9b) 

Five  adaptive  and  uniform  mesh  solutions  of  Eq.  (8)  were  calculated  for  0  <  r  S  3.2 
and  our  findings  are  summarized  in  Table  1.  Solutions  3  and  4,  with  refinement,  were  cal¬ 
culated  using  an  error  tolerance  of  0.0002  and  a  maximum  of  two  levels  of  refinement. 
The  tolerance  and  maximum  level  of  refinement  were  selected  so  that  the  high-error  region 
under  the  cone  would  maintain  approximately  the  same  mesh  spacing  as  the  uniform  mesh 
used  to  obtain  Solution  5,  The  grids  that  were  used  to  obtain  Solution  4  are  shown  in 
Figure  5  at  f  =  0.56,  1.68,  2.24,  and  3.2.  A  new  base  mesh  was  introduced  at  f  =  2.82. 
The  meshes  that  were  used  to  obtain  Solutions  2,  3,  and  4  at  r  =  3.2  are  shown  in  Figure 
6.  Finally,  surface  and  contour  plots  of  Solutions  1,  2,  and  3  and  of  Solutions  4  and  5  at 
r  =  3.2  are  shown  in  Figures  7  and  8,  respectively. 

Solution  1  bears  no  resemblance  to  the  exact  solution  and  demonstrates  the  devastat¬ 
ing  effects  of  large  dissipative  and  dispersive  errors.  Solution  2,  with  mesh  moving  only, 
provides  a  dramatic  improvement  in  the  results  for  approximately  one-half  the  cost  of 
using  both  mesh  motion  and  refinement.  Solution  5  took  more  than  three-times  longer  to 
calculate  than  Solution  4  for  approximately  the  same  accuracy;  thus,  demonstrating  the 
efficiency  of  the  refinement  process.  The  subgrids  for  the  refined  Solutions  3  and  4  are 
concentrated  in  the  region  of  the  cone  and  are  aligned  with  its  principal  axes  as  it  rotates. 


Strategy 

Base 

Mesh 

Base  Time 
Step 

Ik  111 

IklL 

CPU  Time 
(sec.) 

Stationary 
uniform  mesh 

14x14 

0.056 

0.2560 

0.78 

46 

Moving  mesh 

32x32 

0.026 

0.0301 

0.20 

458 

Stationary  mesh 
with  refinement 

14x14 

0.056 

0.0832 

0.48 

852 

Moving  mesh 
with  refinement 

14x14 

0.056 

0.0249 

0.18 

904 

Stationary 
uniform  mesh 

56x56 

0.014 

0.0759 

0.48 

2647 

Table  1.  Errors  at  r  =  ?.2  and  computational  costs  for  five  solutions  of  Example 
1. 

Dissipative  and  dispersive  errors  cause  a  "wake"  of  spurious  oscillatory  information  to  fol¬ 
low  the  moving  cone  (cf.  Figures  7  and  8).  Some  mesh  refinement  is  performed  in  the 
wake  region  and  this  greatly  reduces  the  magnitude  of  the  oscillations. 

Example  2.  Consider  the  uncoupled  linear  initial- boundary  value  problem 


Ui  +Ui  =0,  Ml  -  “2.  =  ^  >  0.  (xo’)eQ, 


1  -  16((x-l/2)^  -t-  1.5y2),  if  {x-mf  +  1.5y^  ^  1/16 
U\ix,y,0)-  otherwise. 


(x,y)€Q^9n,  (10b) 


J  1  -  16((x+1/2)2  +  1.5y2)^  if  (X -(-1/2)2  +  1.5y2  < 

“2<^’>’0)=  lo,  otherwise. 


(x,y)  e  (10c) 
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=  U2(x,y,t)  =  0,  r  >  0,  (j:,y)  e  dQ,  (lOd) 

and  Q  :=  { (jc .y)  I  -1  ^  j:  ^  1,  -0.6  i^y  S  0.6 ). 

The  solution  of  this  problem  consists  of  two  moving  cones  that  collide  and  pass 
through  each  other.  We  selected  it  in  order  to  determine  how  the  various  adaptive  stra¬ 
tegies  could  cope  with  interacting  phenomena. 

One  uniform  mesh  and  three  adaptive  solutions  of  Eq,  (10)  were  calculated  for 
0  <  r  ^  1.2  and  our  findings  are  summarized  in  Table  2.  The  solutions  involving 
refinement  were  computed  with  a  tolerance  of  0.0038.  All  solutions  were  designed  to 
have  approximately  the  same  accuracy.  The  grids  that  were  used  to  obtain  Solution  4  are 
shown  in  Figure  9  at  r  =  0,  0.23,  0.46,  0.92,  and  1.2. 

The  results  of  Table  2  demonstrate  the  efficiency  of  the  mesh  moving  strategy  on  this 
example.  Solution  2  with  mesh  moving  was  slightly  more  accurate  than  Solution  1 
obtained  on  a  uniform  mesh,  and  it  required  less  than  one-half  of  the  computation  time. 
Solution  3  with  refinement  on  a  stationary  mesh  shows  only  a  modest  improvement  over 
Solution  1;  however,  the  combination  of  mesh  moving  and  refinement  computed  in  Solu¬ 
tion  4  again  shows  a  significant  gain  in  efficiency.  We  suspect  that  the  high  accuracy 
achieved  by  mesh  moving  on  this  example  is  due  to  the  reduction  in  dispersive  errors  that 
results  when  the  mesh  follows  the  cones  with  approximately  the  correct  velocity. 

Example  3.  Consider  the  Euler  equations  for  a  perfect  inviscid  compressible  fluid 


,  ^ 


-28- 


vS 

r 

i 

I 

r 

ishn 

<^kri 


u,  +f,(u)  +  g„(u)  =  0. 


(11a) 


where 


,  f(u)  = 


pu^+p 

puv 

u(e+p) 


g(u)  = 


pv^+p 

vie+p) 


(llb,c,d) 


Here,  u  and  v  are  the  velocity  components  of  the  fluid  in  the  x  and  y  directions,  p  is  the 
fluid  density,  e  is  the  total  energy  of  the  fluid  per  unit  volume,  and  p  is  the  fluid  pressure. 
For  an  ideal  gas 


P  -  (Y- 1)[«  -  P(M^  +  v2)/2],  ( 

where  y  is  the  ratio  of  the  specific  heat  at  constant  pressure  to  that  at  constant  volume. 


We  solve  a  problem  where  a  Mach  10  shock  in  air  (Y=  1.4)  moves  down  a  channel 
containing  a'  wedge  with  a  half-angle  of  thirty  degrees.  This  problem  was  used  by  Wood¬ 
ward  and  Collela  [32]  to  compare  several  finite  difference  schemes  on  uniform  grids. 
Like  them,  we  orient  a  rectangular  computational  domain,  -0.3  <3.4,  0  <1,  so 

that  the  top  edge  of  the  wedge  is  on  the  bottom  of  the  domain  in  the  interval  y  =  0, 
1/6  <  X  <  3.4.  Thus,  in  the  computational  domain  it  appears  like  a  Mach  10  shock  is 
impinging  on  a  flat  plate  at  an  angle  of  sixty  degrees.  The  initial  conditions  that  are 
appropriate  for  this  situation  are 


p  =  8.0,  p  =  116.5,  e  =563.5,  u  =4.125V3,  v  =-4.125, 

if  y  <  V3(jc-l/6), 


m 


Along  the  left  boundary  (x  =  -0.3)  and  the  bottom  boundary  to  the  left  of  the  wedge 
(y  =  0,  -0.3  <x  ^  1/6),  we  prescribe  Dirichlet  boundary  conditions  according  to  Eq.  (12); 
along  the  top  boundary  (y  =  1),  values  are  prescribed  that  describe  the  exact  motion  of  an 
undisturbed  Mach  10  shock;  along  the  right  boundary  (x  =  3.4),  all  normal  derivatives  are 
set  to  zero;  and  along  the  wedge  (y  =  0,  1/6  £  x  ^  3.4),  reflecting  boundary  conditions  are 
used. 

The  solution  of  this  problem  is  a  complete  self-similar  structure  called  a  double-Mach 
reflection  that  was  described  in  Ben- Dor  and  Glass  [12,  13].  Two  reflected  Mach  shocks 
form  with  their  associated  Mach  stems  and  contact  discontinuities.  The  geometry  of  these 
structures  is  very  fine  and  is  primarily  confined  to  a  small  region  that  moves  along  the 
wedge  with  the  incident  shock.  One  of  the  two  contact  discontinuities  is  so  weak  that  it  is 
usually  not  noticed  in  computations. 

The  MacCormack  finite  difference  scheme  needs  artificial  viscosity  to  "capture" 
shocks  without  excessive  oscillations.  We  used  a  model  developed  by  Davis  [20]  which  is 
total  variation  diminishing  in  one-space  dimension. 

Five  solutions  of  this  problem  were  calculated  for  0  <  r  <  1.9  as  indicated  in  Table 
3.  Refinement  was  restricted  to  a  maximum  of  two  levels  and  a  tolerance  of  0.6  in  the 
maximum  norm  was  prescribed.  A  pointwise  error  indicator  based  on  the  assumption  of 
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Strategy 

Base 

Mesh 

Ik  111 

IkIL 

Stationary 
uniform  mesh 

64x34 

0.066 

0.26 

Moving  mesh 

44x20 

0.056 

0.18 

Stationary  mesh 
with  refinement 

44x20 

0.055 

0.23 

Moving  mesh 
with  refinement 

44x20 

0.039 

0.16 

CPU  Time 
(sec.) 


Table  2.  Errors  at  r  =  1.2  and  computational  costs  for  four  solutions  of  Example 

2. 


smooth  solutions,  like  the  present  one,  is  not  appropriate  for  problems  having  discontinui¬ 
ties.  Without  restricting  the  maximum  level  of  refinement,  we  could  refine  indefinitely  in 
the  vicinity  of  a  discontinuity. 

Solutions  2  through  5  were  intended  to  be  of  comparable  accuracy  and  we  shall 
attempt  to  appraise  the  computational  cost  of  each  adaptive  strategy.  The  maximum 
number  of  nodes  that  was  introduced  in  any  base  time  step  and  the  total  CPU  time  are 
presented  as  measures  of  computational  complexity  in  Table  3.  Contours  of  the  density  at 
t  =  0.19  are  shown  for  all  five  solutions  in  Figure  10  and  the  grids  that  were  generated  for 
Solution  4  at  r  =  0.038,  0.076,  0.114,  0.152,  and  0.19  are  shown  in  Figure  11. 

As  in  the  previous  two  examples,  the  mesh  moving  strategy  of  Solution  2  does  a 


Strategy 

Base 

Mesh 

Max.  No. 

Nodes 

CPU  Time 
(sec.) 

Stationary 
uniform  mesh 

63x29 

1827 

2130 

Moving  mesh 

63x29 

1827 

2220 

Stationary  mesh 
with  refinement 

29x11 

2782 

3254 

Moving  mesh 
with  refinement 

29x11 

3540 

3725 

Stationary 
uniform  mesh 

120x40 

4800 

6861 

Table  3.  Maximum  number  of  nodes  in  any  base  time  step  and  computational 
costs  for  five  solutions  of  Example  3. 


great  deal  to  improve  the  results  of  the  static  Solution  1  for  approximately  a  five-percent 
increase  in  computational  cost.  Comparing  the  top  two  contours  oi  Figure  10,  we  see  that 
the  resolution  of  the  incident  and  reflected  shocks  is  much  finer  with  Solution  2  than  with 
Solution  1.  Additional  detail  of  the  structures  in  the  Mach  stem  region  and  of  the  contact 
discontinuities  is  present  in  Solution  2,  but  not  in  the  nonadaptive  Solution  1.  Finally, 
Solutions  1  and  5  display  more  oscillatory  behavior  behind  the  incident  shock  near  the 
upper  boundary.  This  is  undoubtedly  due  to  our  maintaining  a  discontinuity  where  the 
shock  intersects  the  upper  boundary. 

The  use  of  refinement  on  a  stationary  mesh  again  does  not  give  the  dramatic 
improvement  obtained  by  mesh  moving  (cf.  the  second  and  third  contours  of  Figure  10). 


N- ////A* 


Initally  the  fine  meshes  were  following  the  incident  and  reflected  shock  structures  and 


better  results  were  obtained;  however,  by  t  =0.19  refinement  is  being  performed  over 
much  of  the  domain  and  two  levels  of  refinement  are  not  sufficient  for  adequate  resolution 
(cf.  Amey  and  Flaherty  [6]).  The  combination  of  mesh  motion  and  refinement  depicted  by 
Solution  4  in  Figure  10  provides  a  marked  improvement  in  resolution.  The  sequence  of 
meshes  shown  in  Figure  11  shows  that  the  coarse  mesh  is  able  to  follow  the  differing 
dynamic  structures  and  that  refinement  is  only  performed  in  the  vicinity  of  discontinuities. 
Initially,  only  one  rectangular  cluster  is  needed  to  follow  the  incident  shock  (cf.  Amey 
and  Flaherty  [5]).  As  time  progresses,  two  clusters  are  created  in  order  to  follow  the 
incident  and  reflected  shocks  (cf.  the  upper  three  meshes  of  Figure  11).  A  third  cluster  is 
created  as  time  increases  funher  in  order  to  follow  the  evolving  activity  in  the  region  of 
the  Mach  stem  (cf.  the  lower  two  meshes  of  Figure  11). 

Severe  distortion  of  the  mesh  in  the  reflected  shock  region  caused  a  static  mesh 
regeneration  to  occur  for  Solution  4  at  r  =  0.162.  The  base  meshes  before  and  after  the 
static  regeneration  are  shown  in  Figure  12.  Thus,  Solution  4  demonstrates  all  of  the  capa¬ 
bilities  of  our  adaptive  procedure.  The  results  presented  in  Table  3  and  Figure  10  also 
show  that  Solution  4  provided  greater  resolution  than  the  uniform  mesh  Solution  5  for 
approximately  one-half  of  the  cost.  Solution  4  also  shows  many  of  the  same  characteris¬ 
tics  as  the  solution  computed  by  Woodward  and  Collela  [32]  using  MacCormack’s  method 
on  a  240  X  120  uniform  grid.  We  were  unable  to  compute  a  solution  on  such  a  fine  mesh 
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due  to  virtual  memory  limitations  on  our  computer,  however,  we  estimate  that  it  would 
have  used  14,400  nodes  and  40,000  CPU  seconds. 


Figure  12.  Base  grids  before  (top)  and  after  (bottom)  the  static  mesh  regenera¬ 
tion  that  was  performed  for  Solution  4  of  Example  3  at  r  =  0.162. 

The  results  presented  for  this  problem  demonstrate  the  power  and  efficiency  of  our 
adaptive  techniques;  however,  we  would  have  preferred  to  allow  more  than  two  levels  of 
refinement  and  a  finer  base  mesh.  These  calculations  would  have  produced  better  resolu¬ 
tion  of  the  discontinuities  and  other  fine-scale  structures  that  further  demonstrate  the 


computational  advantages  of  adaptive  methods  relative  to  uniform  mesh  techniques.  As 


noted,  restrictions  of  our  computing  environment  prevented  us  from  doing  this  in  a  reason¬ 


able  manner.  We  hope  to  perform  these  calculations  in  the  future  using  a  larger  comput¬ 


ing  system. 


IV.  DISCUSSION  OF  RESULTS  AND  CONCLUSIONS. 


We  have  described  an  adaptive  procedure  for  solving  systems  of  time-dependent  par¬ 


tial  differential  equations  in  two-space  dimensions  that  combines  existing  mesh  moving  [5] 


and  local  refinement  [6]  techniques.  The  algorithm  also  contains  procedures  for  initial 


mesh  generation  and  static  mesh  regeneration.  It  can  be  used  with  a  wide  variety  of  finite 


difference  or  finite  element  schemes  and  error  indicators. 


We  obtained  computational  results  for  hyperbolic  systems  of  conservation  laws  by 


using  our  adaptive  methods  with  a  MacConnack  finite  difference  scheme  and  using 


Richardson’s  extrapolation  to  furnish  local  error  indicators.  Our  computational  results  on 


three  examples  indicate  that  mesh  moving  can  significantly  reduce  errors  for  approxi¬ 


mately  a  ten-percent  increase  in  cost  relative  to  computations  performed  on  stationary  uni¬ 


form  meshes.  The  use  of  local  refinement  without  mesh  moving  provided  increased 


efficiency  relative  to  uniform-mesh  calculations,  although  not  as  dramatic  as  those  found 


using  mesh  moving.  The  combination  of  mesh  moving  and  local  refinement  provided  reli¬ 


able  results  while  costing  significantly  less  than  stationary-mesh  calculations.  Thus,  the 


mmmm 


overhead  associated  with  the  dynamic  data  structures  is  less  than  the  time  to  calculate  a 
comparable  solution  on  a  uniform  mesh. 

The  results  of  Section  HI  and  others  (cf.  Amey  and  Flaherty  [5,  6])  indicate  that  our 
mesh  moving  procedures  perform  better  alone  than  with  refinement.  This  is  because  the 
projection  of  fine-mesh  solutions  onto  coarser  meshes  reduces  the  errors  at  base  mesh 
nodes,  and  mesh  motion  based  on  controlling  small  or  zero  local  discretization  errors 
either  fails  or  results  in  no  movement.  Erratic  mesh  motion  can  also  occur  with  some 
techniques  when  movement  indicators  are  small.  This  topic  is  discussed  in  Coyle  et  al. 
[19]  and  a  possible  remedy  for  one-dimensional  problems  is  suggested  in  Adjerid  and 
Flaherty  [2],  Further  experimentation  and  analysis  are  being  performed  in  order  to  deter¬ 
mine  the  best  way  to  combine  mesh  moving  and  refinement. 

There  are  several  other  ways  to  improve  the  efficiency,  reliability,  and  robusmess  of 
our  adaptive  methods.  The  present  Richardson’s  extrapolation-based  error  indicator  is 
expensive  and  we  are  seeking  ways  of  replacing  it  by  techniques  using  p-refinement. 
Such  methods  have  been  shown  [1,  2,  3,  8,  10,  16,  22]  to  have  an  excellent  cost  perfor¬ 
mance  ratio  when  used  in  conjunction  with  finite  element  methods.  An  appropriate  error 
indicator  or  estimator  can  be  used  to  control  a  differential  refinement  algorithm,  where 
different  refinement  factors  (i.e.,  other  than  binary)  are  used  in  different  high-error  clus¬ 
ters.  If  the  error  indicator  is  capable  of  providing  separate  estimates  of  the  spatial  and 
temporal  errors,  as  the  present  one  does,  then  different  refinement  factors  can  also  be  used 


in  space  and  rime.  We  also  hope  to  demonstrate  the  flexibility  of  our  refinement  pro¬ 


cedure  by  using  it  with  a  finite  difference  or  finite  element  scheme  for  parabolic  problems. 

The  greater  reliability  and  efficiency  of  adaptive  techniques  will  be  most  beneficial  in 
three  dimensions.  These  techniques  must  be  able  to  take  advantage  of  the  latest  advances 
in  vector  and  parallel  computing  hardware.  The  tree  is  a  highly  parallel  structure  and  we 
have  been  developing  solution  procedures  that  exploit  this  in  a  variety  of  parallel  comput- 
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